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' Abstract 

. In this paper we make a detailed numerical comparison between three algorithms 

for the computation of the full Lyapunov spectrum as well as the associated eigen- 



o 



' vectors of general dynamical systems. They are : (a) the standard method, (b) 

On , a differential formulation of the standard method, and (c) a new algorithm which 

does not require rescaling and reorthogonalization. We also bring out the relations 
among these methods. Moreover, we give a simplified formulation of the new algo- 
■ rithm when the dimensionality of the system is 4. We find that there is reasonable 

Q , agreement among the Lyapunov spectra obtained using the three algorithms in most 

I cases. However the standard method seems to be the most efficient followed by the 

1^ ■ new method and the differential version of the standard method (in that order), 

as far as the CPU time for the computation of the Lyapunov spectra is concerned. 
The new method is hardly suitable for finding the eigenvectors, whereas the other 
^ ' procedures give nearly identical numerical results. 

PACS numbers: 05.45.+b, 02.20. Qs 

1 Introduction 

Extreme sensitivity to initial conditions is the commonly accepted defining property of 
chaos in nonlinear systems. Lyapunov exponents which determine the exponential rates 
at which nearby trajectories diverge on an average, are the quantitative characteristics of 
a chaotic orbit. A dynamical system of dimension n has n Lyapunov exponents and n 
principal directions or eigenvectors, corresponding to a set of nearby trajectories [1]. One 
of the standard and popular methods to compute the Lyapunov spectrum of a dynamical 
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system involves a Gram-Schmidt Reorthonormalizaton (GSR) of the 'tangent vectors' 
[2]. A differential version of this method has been formulated which corresponds to a 
continuous GSR of the tangent vectors [3]. A modification of this method with the 
introduction of a stability parameter makes it dynamically stable, applicable to systems 
with degenerate spectra, and reliable for computations [4]. Recently, a new algorithm 
for the computation of Lyapunov exponents has been proposed, which has been claimed 
to be vahd even for evaluating partial Lyapunov spectra [5]. This is based on the 'QR' 
method for the decomposition of the tangent map, (where Q is an orthogonal matrix and 
R is an upper triangular matrix) which has been studied by several authors [6] . It utilizes 
representations of orthogonal matrices applied to the tangent map, and does not require 
the GSR procedure. It has also been claimed that it has several advantages over the 
existing methods, as it involves a minimum number of equations. In this paper we have 
made a detailed comparison of the three algorithms as regards accuracy and efficiency, 
by computing the full Lyapunov spectra of some typical nonlinear systems with 2,3 and 4 
variables. We also compare the performance of the standard method with its differential 
version, in computing the Lyapunov eigenvectors. 

In section 2, we outline the three methods with necessary details. We bring out the 
relation between the differential version of the standard method and the new procedure, 
by deriving the differential equations of the latter from those of the former. It is difficult 
to use the new method with a standard representation of orthogonal matrices when the 
number of dimensions of the system is greater than 3. In section 3, we give a convenient 
representation for them for n = 4, by making use of the well-known fact that 5'0(4) ~ 
5*0(3) X S0{3) [7]. This simplifies the calculations considerably. In section 4, we make 
a comparative study of the three algorithms for the computation of Lyapunov spectra 
by taking up some typical 2, 3 and 4 dimensional systems. We have considered both 
dissipative and Hamiltonian systems of some physical interest, for comparison. In section 
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5, we compare the computation of the Lyapunov eigenvectors (which are local properties), 
using these algorithms. In section 6, we make a few concluding remarks. 

2 Computation of Lyapunov exponents 

Consider an n-dimensional continuous-time dynamical system: 

f = F(Z,t). (1) 

where Z and F are n-dimensional vector fields. To determine the n Lyapunov exponents 
of the system, corresponding to some initial condition Z(0), we have to find the long term 
evolution of the axes of an infinitesimal sphere of states around Z(0). For this, consider 
the tangent map given by the set of equations, 

where J is the n x n Jacobian matrix with 



J--^ (3) 



A solution of equation (2) can be formally written as 

5Z(t) = M(Z(t),t)5Z(0), (4) 
where M(Z(t), t) is the tangent map whose evolution equation is easily seen to be 

— =J.M. (5) 

In the following, we give a brief description of the procedures for computing the n 
Lyapunov exponents of the system using (a) the standard method, (6) the differential 
version of the standard method and (c) the new method based on the 'QR' decomposition 
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of M, which dispenses with the tangent vectors 5Z, and in a sense, computes the exponents 
directly. 

(a) The Standard method 

Let Ai, A2, . . . , A„ be the n Lyapunov exponents of the system in a decreasing sequence, 
Ai > A2 . . . > An- In the standard method [2] one first chooses n orthogonal tangent 
vectors as initial conditions for eq.(2). The standard choice is ei(0) = (1, 0, 0, . . .); 62(0) = 
(0, 1, 0, 0, . . .), etc.. Eq.(2) is then solved upto time r for each of the initial conditions 
yielding vectors vi(r), V2(t), . . . Vn('r). These vectors are orthonormahzed using a Gram- 
Schimdt Reorthonormalization (GSR) procedure to yield: 



e.(r) = -^-(-^-^^(-))^^(-) (6) 
||V2- (V2,ei(r))ei(r)|| 

and so on. The norms in the denominators, denoted by A^i(l), A^2(l), ■ ■ ■ -^n(l)) are 
stored for the computation of Lyapunov exponents. The procedure is repeated for a sub- 
sequent time r of integration using ei(T) as initial conditions for eqn.(2). The resulting 
vectors Vi(2r), are again orthonormahzed using a GSR procedure to yield orthonormal 
tangent vectors ei(2T), i — and the norms Ni{2), N2{2)^ . . . Nn{2). After r it- 

erations, we get the orthonomal set of vectors ei{rT),i — 1,. . . ,n at time t — rr. The 
Lyapunov exponents are 

This is due to the following reason. Since GSR never affects the direction of the first 
vector in a system, this vector tends to seek out the direction in the tangent space, which 
is most rapidly growing and its norm is proportional to e^^* for large t. The second vector 
has its component along the direction of the first vector removed and its norm would be 



proportional to e'*'^* for large t and so on. 

It is to be noted that we have to integrate n{n + 1) coupled equations in this method, 
as there are n equations for the fiducial trajectory in (1) and n copies of the tangent map 
equations in (2). 

(b) Differential version of the standard method 

In this method [3], the orthonormal set of vectors ei{t) are obtained by solving dif- 
ferential equations set up for them, instead of resorting to the GSR at discrete steps. 
Rather, GSR is incorporated in the procedure itself. It can be shown that 

d '-^ 

— ej(t) = Gsi - GiiBi - YliGij + Gji)ej, (8) 
where G = J is the Jacobian matrix introduced in eq.(2) and 

G,j^{e,{t),J{Z{t))e,{t)), (9) 

that is, Gij are the matrix elements of the Jacobian in the basis ei{t). Now let 6^(0) 
evolve to ei{t). 

e,(t) = M(Z(t),^)e,(0), (10) 
In fact, ei{t) is the orthonormalized set corresponding to ei{t) i = 1, .. .n. Define 

d,j = {e,{t),e^{t)). (11) 
The GSR procedure ensures that dij is a lower triangular matrix: 

dij^O, i<j. (12) 

It can be shown that 
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n, 



(13) 



and that, 




(14) 



for large t. That is, 



Xi = hm - log da. 



(15) 



The Lyapunov exponents are computed by solving the coupled equations (1), (8) 
and (13) in this method. As there are equations for the n components each of the 
orthonormal vectors ei{t) in (8), n equations for da in (11), apart from the n equations 
for the fiducial trajectory in (1), we have to integrate n{n + 2) coupled equations in this 
method. 

In practice, this procedure is not numerically 'stable', as the set ei{t) may not remain 
orthonormal under the time evolution. In particular Ajj = (ej(t), ejit)) — Sij, 1 < i, j < n 
may not all vanish. Moreover, the method is not applicable to systems with degenerate 
exponents. These are remedied by a modification of the method, using a stability pa- 
rameter /3 [4]. We replace Gu by Gu + (3{{ei,ei) — 1) and Gij by Gij + I3{ei,ej),i ^ j 
in equations (8) and (13). Though it has been shown that the method is strongly stable 
when (3 > — A„, where A„ is the lowest exponent, it is found in certain problems, that 
(3 has to be significantly larger than — A„ in practice. Moreover, it may be pointed out 
that this method requires prior knowledge of the lowest Lyapunov exponent for the 
computation of the complete spectrum Aj. If an arbitrarily high value is assigned to 
one ends up with an arithmetic overflow problem during computations. 



(c) The New method based on a 'QR' decomposition of M 
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The new algorithm [5] is based on a 'QR' decomposition of M, where Q is an orthogonal 
matrix and R is an upper triangular matrix. This results in a set of coupled differential 
equations for the Lyapunov exponents along with the various angles parametrising the 
orthogonal matrices. In this subsection we derive these equations from the differential 
version of the standard method considered in the previous subsection . 

Consider the tangent map matrix M. From eq.(lO), 

= (e,(0),Me,(0)) = (e,(0), e,(t)) (16) 
As ej{t) form an orthonormal set of vectors, we have from eq.(ll), 

ej{t)^J2^k{t)djk. (17) 

Hence, 

Mij^^{ei{Q),e„it))djk. (18) 

k 

Define the matrices Q and R by 

Qij = {em,ej{t))^{ej{t))„ (19) 
and 

Rij = dji. (20) 

Hence, 



M^QR. (21) 

Clearly the columns of Q are the orthonormal vectors ej{t), and Q is an orthogonal 
matrix. As (i is a lower triangular matrix, R is an upper triangular matrix. 
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Now Gij and Jij are the matrix elements of the Jacobian in the orthonormal bases 
ei{t) and ei(0) respectively and related by a rotation transformation represented by Q. 

Introducing complete sets of states at the appropriate places, we have 



= ^(e,(0,efc(0)) (efc(O), JeKO))(eKO),e,(0) 

k,l 

— y^jQikJklQli — {QJQ)ii- (22) 

k,l 

Taking the scalar product of eq.(8) with 6^(0) and making appropriate changes of 
indices, we have 



d d 

J^Qjk = ^(ei(0),efc(t)) 

fc-i 

= (e,(0), Jefe(t)) - Gkk{eMMt)) - Y^i^^i + Gik){ej{0),ei{t)) 

1=1 

k-l 

= {ej{0),Jek{t))-GkkQjk-J2iGki + Gik)Qji. (23) 

1=1 

As all the quantities are real, 

Qij = Qji = (e,(0),ei(t)) = {ei{t),ej{0)). (24) 
Multiplying eq.(23) by Qij on the right and using the fact that 



%(e,(0),Je,(t)) = E(e.(^),e,(0))(e,(0),Je,(i)) 

j 

= {ei{t),3ek{t)) = Gik (25) 

we find 
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= Gik — Gkk X] QijQjk — X] ^i^kl + Glk)QijQjl 
3 3 1=1 

k-1 

— Gik — GkkSik — ^^{Gki + Gik)6ii, (26) 
1=1 



as Q is an orthonormal matrix. 



Again, Q-^Q is an antisymmetric matrix as Q is orthogonal and it is sufficient to 
consider i > k. In tfiis case, tfie last term vanishes and we obtain. 



(Qj^Q)ik = = (QJQhk, i > k. (27) 

Q being an orthogonal matrix is characterised by angles and we obtain differ- 

ential equations for these angles. Prom eqs.(13) and (14), the differential equations for 
the Lyapunov exponents are 

l(X,t) = Gu = {QJQh. (28) 

In this method, we have essentially traded the orthonormal vectors ej(t) for the orthog- 
onal matrix Q parametrized by the angles. We have to solve the coupled eqs.(l), 
(27) and (28) in this procedure to obtain the Lyapunov exponents. We have to integrate 
n + + n = ^^^^^^ coupled equations in this method. 

3 A convenient representation for Q and simplifica- 
tion of QQ for n = 4 

In [5] , the explicit representation of the orthogonal matrix Q used is the one in which it 
is represented as a product of '^^^^^^^ orthogonal matrices, each of which corresponds to a 
simple rotation in the i — j*^ plane {i < j) . Thus Q 



9 



Q ^ Q(12)Q(13)Q(14) _ _ _ Q{ln)Q{2Z) _ _ _ Q{n-2,n-l) Q{n-l,n) 

where 

= cos 6ij if k = I = iorj; (29) 

= sin 6ij if k = i,l = j; 

— — sin % if k — j,l — i; 

— otherwise. 

In terms of the group generators, O^*-'^ can be written as 

Q{ij) ^ Q0ij{tij)^ (30) 

where the generator tij is represented by 

{'tij)kl = ^ik^jl — ^il^jk, (31) 

The generators satisfy the commutation relations, 



[^ijj^mn] — ^intjm ~l~ ^jm^in ^imtjn ^jn^im- (32) 

The above representation for Q is conceptually simple and works very well for n = 
2 and 3 [5]. However, for n > 3, it is hardly suitable for practical computations of 
Lyapunov exponents. This is because the expressions for QQ and QJQ are very lengthy 
and unmanageable even for n = 4. 

In the present work, we employ a representation for Q, which simplifies the calculations 
and numerical computations for n = 4. This is based on the well-known fact that SO {4) ~ 
S0{3) X 50(3) [7]. Prom the generators tij we construct the following combinations: 

Ml = 1(^23 + tu), Ni = l(t23 - tu) 

M2 = |(% + i24), iV2 = !(% - ^24) (33) 

Ms = i(ti2 + t34), ATg = i(tl2-^34)■ 
Then it is easily verified that Mj and generate two mutually commuting SO {3) 

algebras: 
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(34) 



We write Q as 



Q — QiiQi 



(35) 



where 



(36) 



and 



(37) 



Using 



e^Fe-^ = y + [X, y] + -[X, [X, Y]] + ..., 



(38) 



for any matrices X and F and the commutation relations in equation (32), it can be easily 
verified that 

QQ =QiQi + QnQii 

= [01 + 03 sin 62] Ml + [62 cos 9i + 63 sin 0i cos 62] M2 (39) 

+[62 sin 02 + 63 cos 6'i cos 62^3 + [64 + 9^ sin 6*5] A^i+ 

+[^5 cos 64 — 6*6 sin 64 cos ^5] A^2 + [^5 sin 9^ + 9% cos ^4 cos 9^]N3. 

The exphcit form of the matrices MjandA^j can be found using equations (31) and (33) 
and are written in terms of 2 x 2 blocks as given below: 



-(7i 



-1(72 



(7l 




-^c^2 




Mo 



(73 



-0-3 



M3 



l(J2 



^c^2 








-1(72 



(40) 
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Here I is the 2x2 identity matrix and ui, (J2 and cxs are the Pauh matrices: 



"10" 




" 


1 " 




" -i " 




" 1 " 


1 


, 0-1 = 


1 





, (72 = 


i 


, 0-3 = 


-1 



Then we find that 



QQ 



I 



-fi{0,9] 



MO, 9) 



-f2i9,9] 
-U{9,9] 


h{9,9) 



-h{9,9) 
-M9,9) 
-h{9,9) 




(41) 



(42) 



where 



/i = —\{92 sin 9i + 9^ cos 9i cos 92 + ^5 sin 9^ + cos ^4 cos 6*5), 

/2 = \ {92 cos 6'i — 9^ sin ^1 cos ^2 + ^5 cos 6*4 — 9q sin 6*4 cos 6^5) , 

/a =-|(^i + e3sin^2-^4-4sin^5), 

/4 =- 1(^1 + ^3 sin ^2 + ^4 + ^6 sin ^5), 

/5 = |(— ^2 COS ^1 + 93 sin ^1 cos ^2 + ^5 cos ^4 — 6*6 sin ^4 cos ^5), 

/e = — I (^2 sin 9i + 9^ cos ^1 cos ^2 — ^5 sin ^4 — cos ^4 cos ^5) . 

Using equation(12), we find that the equations for 9i, spht neatly into two sets: 



/ -1 -Sin^2 \ ( ( G'32 + G41 \ 

-sin^i -cos^iCOS^2 ^2 = G'2i + G'43 , (44) 

»^ 008 6*1 -sin9i cos 92 J \ 9^ J \ G31-G42 J 

and 





-1 


— sin 9^ \ 


( \ 


( G32 


— G41 




— sin 6*4 


—cos94 cos ^5 


k = 


G21 


— G43 


V 


cos ^4 


—sin9i cos 6*5 / 


\ 9, 1 


[ G31 


+ G42 



We also have 

j^{X,t)^G,„i^l,...A, (46) 

fi:om eq.(28). Hence, to find the Lyapunov exponents of a dynamical system with 4 
variables, we have to solve the evolution equations for the system given by eq.(l) and the 
tangent map equations given by equations (44), (45) and (46) after finding G = QJQ. 

Any 4x4 matrix J can be written as : 



16 

J = J2aiXi, (47) 

i=l 

where the 16 matrices are defined in terms of 2 x 2 blocks as 



12 



X 



13 



/ 
I 

ai 

(7i 

/ 
/ 



(7i 
(71 



-^10 — 



X 



14 



-/ 
/ 

-ai 

(7i 

/ 
-/ 



(7i 
-(7l 



15 



(T3 

(73 

i(72 

i(72 

" (73 " 
(T3 



— i(72 
-i(72 



X 



12 



16 



-(73 
(73 

i(72 
— i(72 

" (73 " 
-^3 



— i(72 
i(72 



(48) 

It is easy to find commutators [Xi, Mj] and [Xi,7Vj] from eqs.(40) and (48). Then, 
using eqs. (35-38), we can obtain 



G = QJQ. (49) 

4 A comparative study of the three algorithms for 
the computation of Lyapunov spectra 

The standard algorithm involves an explicit GSR for finding the orthonormal set ei{t) 
and the Lyapunov spectrum. The differential version considered in section 3(b) amounts 
to computing the spectrum with continuous GSR. Here explicit GSR is avoided as it is 
incorporated in the method. However, the differential equations for ei{t) in this method 
are nonlinear, as they involve (ej(t), Jej(t)) in the RHS, in contrast to the standard 
method which uses the linearized equations for SZ directly. In the new method, one 
deals directly with the orthogonal matrix relating ei{t) and 6^(0). It uses a minimal 
number of variables and rescaling and reorthogonalization are eliminated. However, in 
this method, the evolution equations for the angles and Lyapunov exponents are highly 
nonlinear involving sines and cosines of the angles. Hence it is not clear 'a priori' which 
method is 'superior' and there is a need to compare the efficiency and accuracy of the three 
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methods. That is the subject matter of the present investigation. Here we consider some 
typical nonhnear systems of physical interest with n — 2,3 and 4. The driven Van der 
Pol oscillator is taken as an example of a two dimensional system, whereas the standard 
Lorenz system is chosen for = 3. For n — 4, we consider the coupled quartic oscillators 
and anisotropic Kepler problem as examples of conservative Hamiltonian systems and the 
Rossler hyperchaos system as an example of a dissipative system. We give the differential 
equations for these dynamical systems in the following. 

1. Driven Van der Pol oscillator {n — 2): 

— ( = ( ^ (50) 

dt\Z2j \ -{d{l- zl)z2- zi + bcosujt, J ' ^ ' 

where b and d are parameters and uj is the driven frequency. In our numerical work we 

have chosen d = —5.0, b = 5.0 and uj = 2.47 as the parameter values. 



2. Lorenz system {n — 3): 

d 



^ zi\ / a{z2 - zi) 

Z2 \ = \ Zi{p - Z3) - Z2 \ . (51) 

\ Z3 J V ^1^2 - f^Zs 

This system is too well-known to require any further discussion. For computations we set 



a = 10.0, p = 28.0 and /3 = |. 

3. Coupled quartic oscillators (n = 4): 

This is a conservative system and the Hamiltonian is given by 

H^^ + ^ + zt + zt + azlzl (52) 

where zi and Z2 are the canonical coordinates, z^ and Z4 are the corresponding mo- 
menta and q; is a parameter. The Hamiltonian in eq.(53) finds applications in high energy 
physics [8] , to mention just one example. The equations of motion are: 

14 



di 



( ^ 




/ 










Zi 










V ^4 J 




[ 


-{4:Z^ + 2azlz2 



\ 



(53) 



This system is known to be integrable for a = 0, 2 and 6 [9]. 
4. Anisotropic Kepler problem (n — A): 
The Hamiltonian of this system is given by: 



(54) 



where 7 is a number. 



The Hamiltonian given above describes the motion of an electron in the Coloumb field 
in an anisotropic crystal, where its effective mass along the x-y plane and z-direction are 
different [10]. 7 = 1 corresponds to the isotropic case and is integrable. When 7 7^ 1, the 
system is non-integrable. Because of the singularity at p = 2; = 0, the Hamiltonian in the 
above form is hardly suitable for numerical integration. For this we choose Zi — p + z 
and Zi = sj p — z as the canonical variables. We can find the corresponding canonical 
momenta ^3 and ^4 in terms of and p^. We also use a re-parametrized time variable r 
defined by (it = (ir(2;f + 2;|). 

The original Hamiltonian with the old variables and energy E corresponds to the 
following Hamiltonian with H' — terms of the new variables [1 1] : 



^' = 2 = \{zl + zl) - E{zl + zl) + (7 - 
The equations of motion resulting from this are: 



(55) 
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( ) 




d 






di 













^3 + (7 -1)^1 ^ 

1^-v (2123-2224) 

Z4 - (7 - 1)^2 ^ 

/ -I ^ {zIzizI+Z2Z3Za{zI-zI)-zIzizI) 
^ ' ~ {z(+zl) 



, 1 \ {zlz2zl-ZlZ3ZA(zl-zl)-zlz2zl) 

We have chosen 7 = 0.61 for computational purposes. 



(56) 



5. Rossler hyperchaos system (n=4): 

This is a dissipative system and an extension of the three dimensional Rossler attractor 
[9]. It is described by the equations: 



di 



( ^1 \ 

\Z4 J 



( -{Z2 + Zz) \ 

zi + az2 + Z4, 
b + Z1Z3 
V CZ4, - dzs J 



(57) 



where a,b,c and d are parameters whose values are taken to be 0.25, 3.0, 0.05 and 0.5 
respectively for our computations. 

In all these cases, the full Lyapunov spectrum is computed using the three methods. 
The time of integration is chosen to ensure reasonable convergence of the Lyapunov expo- 
nents. In most of the cases the time of integration was t — 1, 00, 000 (the exceptions are 
the anisotropic Kepler problem and the Rossler hyperchaos system using the differential 
version of the standard method due to the problem of numerical overflow). For all the 
systems, we have used a variable step-size Runge Kutta routine (RKQC) for integration, 
with an error tolerance, e ~ 10~^ — 10^*. All the computations were performed on a DEC 
Alpha based workstation running OpenVMS. The CPU time taken for each system with 
each of the algorithms was noted. This is the actual time taken by the CPU to accomphsh 
a specific process (independent of the other processes running in the system) . The details 
of the comparison between the two methods are summarized in table 1. 



It may be noticed that all the methods yield essentially the same Lyapunov spectrum. 
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For any autonomous dynamical system, one of the Lyapunov exponents has to be zero 
(corresponding to the difference vector lying along the trajectory itself). For the 
Lorenz system, the Rossler hyperchaos system (both dissipative) and the coupled quartic 
oscillators, this condition is satisfied by all the algorithms. For the anisotropic Kepler 
problem all the methods fail the test. This aspect needs to be studied further. For 
Hamiltonian systems, for every eigenvalue A, there is an eigenvalue —A. This symmetry 
is respected by all the algorithms. For the coupled quartic oscillators, all the exponents 
should be zero corresponding to the integrable case of a = 6. This is indeed satisfied by 
all the algorithms. In Fig.l we give plots of Lyapunov exponents as functions of time, for 
a typical case. Again, there is no significant difference between the three algorithms as 
far as the convergence of the Lyapunov exponents is concerned. It is noteworthy that the 
differential method works well for even systems with degenerate spectra like the coupled 
quartic oscillators. 

On the whole, the standard method seems to have an edge over the new method as 
far as the CPU time for the computation of the Lyapunov sepctrum is concerned. The 
differential version of the standard method generally consumes more CPU time compared 
to the other two methods. For some systems like the anisotropic Kepler problem and the 
Rossler hyperchaos system, there are numerical overflow problems, whatever be the values 
of /3 and the error tolerance e one chooses for this algorithm. In fact, it appears that the 
value of /3 has to be significantly higher than — A„ (indicated by the stability analysis) for 
these systems, for reasonable convergence. 

For the system of coupled quartic oscillators, the CPU time is abnormally high for the 
new method, corresponding to the nonintegrable case of a = 8. This is true both for small 
and large energies. For large energies (~ 25000), since the energy varied by ~ 15 when 
we used the RKQC routine, we also used a symplectic procedure which eliminates secular 
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variations in the energy [13]. With this routine, the CPU times were nearly the same 
for both the methods. However the new method yields poor results for the Lyapunov 
spectrum. For instance corresponding to the initial condition zi — 7.0, Z2 = 7.0, = 5.0 
and Z4 — 4.0, the Lyapunov spectrum computed using the new and the standard methods 
are (1.5506, 0.3254, -0.3261, -1.5499) and (1.5205, 0.0001, -0.0001, -1.5205) respectively. 
The differential version of the standard method led to a numerical overflow problem, 
corresponding to this initial condition. 

In the standard method, after solving for the fiducial trajectory, the equations for the 
tangent flow are linearized equations. In method (b), corresponding to continuous GSR, 
these equations are nonlinear. In the new method, these equations are replaced by the 
equations for the angles determining the principal axes or the bases associated with the 
Lyapunov spectrum and the Lyapunov exponents. These equations involving sines and 
cosines of the angles are highly nonlinear. For dissipative systems this nonlinearity does 
not pose a problem. However in many cases, this nonlinearity renders the differential 
version of the standard method and the new method less efficient and can even lead to 
inaccuracies, in strongly chaotic situations. 

5 Lyapunov eigenvectors 

Earlier we had defined the matrix d,,- as: 




(58) 



Consider the quantities 



'33 



= 0,i<j. 



(59) 
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Define the vectors di as 



di — {dii,d2i,dsi, . . .), 

J2 = {0,d22,d32, ■ ■ ■),etc. (60) 

Let Di, D2, . . . , Dn, be the orthonormal set of vectors obtained form d[s by the Gram- 
Schmidt procedure, starting with di. It can now be shown that Di, D2, . . . , Dn, are 
the eigenvectors of MM or the Lyapunov eigenvectors corresponding to the eigenvalues 
Ai, A2, . . . , A„ [3]. In this section, we consider the compuation and convergence of these 
eigenvectors corresponding to the systems considered in section 4. 

In the standard method, we have to compute ej(t) and ej(t) separately to obtain dij. 
As all the vectors ej{t) tend to ahgn along ei, both dij — ei{t).ej{t) and djj — ej{t).ej{t) 
would tend to zero for j > 1. As dij is the ratio of d^j amd djj, it would be difficult 
to compute them for large t, in this method. Even then, the procedure seems to give 
reasonable results for all the systems, we have considered. 

In the differential version of the standard method, it has been shown that dij satisfy 
the following differential equations [3] : 

d - * dkk 

T+^ij^ -j—dik{Gjk + Gkj),i>j. (61) 

k=j+i %i 

So the eigenvectors are obtained by direct integration of these equations. This pro- 
cedure does not pose any problem as we do not come across division by small numbers 
here. Indeed we find that the eigenvectors converge much more rapidly than the Lyapunov 
exponents in all the cases, as anticipated by Goldhirsch, et al.[3]. 

In the new method, the orthonormal vectors ei{t) are just the columns of the orthog- 
onal matrix Q. However, it is not straightforward to compute ei{t) in this method. So 

19 



we do not consider this method further here. 

We summarise the results for the Lyapunov eigenvectors in table 2 for the same sys- 
tems with the same parameters and initial conditions as in table 1. As remarked earlier, 
the vectors converge sufficiently fast and the two methods yield essentially identical re- 
sults. Now for a Hamiltonian system, the tangent map matrix M satisfies the 'sympletic 
condition': 



MSM = S, (62) 
with 



S = 







(63) 



-/ ; 

where and / are (| x |) null matrix and identity matrix, respectively [11]. It can 
be shown that if D is an eigenvector corresponding to eigenvalue A, then the eigenvector 
corresponding to the eigenvalue —A is SD [1]. This symmetry is very evident in our 
numerical values of the eigenvectors in the case of coupled quartic oscillators, but is 
satisfied only approximately in the case of the highly nonlinear anisotropic Kepler problem. 
It is to be noted that the eigenvectors are dependent upon the initial conditions and are 
only 'local' properties. 

6 Conclusions 

In a recently proposed new method [5] the Lyapunov exponents, are computed directly, 
so to say, by utilizing representations of orthogonal matrices, applied to the tangent map. 
In this paper, we have established the connection between this method and a 'differential 
formulation' of the standard procedure to compute the Lyapunov spectra. We have also 
used the standard decomposition 50(4) ~ S0{3) x 5*0 (3) to simphfy the calculations 
for n = 4, which are otherwise very involved. It has been claimed that the new method 
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has several advantages over the existing methods as it does not require renormahzation 
or reorthogonahzation and requires lesser number of equations. This led us to make a 
detailed comparison of the new method with the standard method as well as its differential 
version , as regards accuracy and efficiency, by computing the full Lyapunov spectra of 
some typical nonlinear systems with 2, 3 and 4 variables. There is reasonable agreement 
among the three procedures as far as the values of the Lyapunov exponents are concerned. 
However, the standard method seems to score over the other two, as far as efficiency (as 
indicated by the CPU time for a process) is concerned, especially in certain strongly 
chaotic situations, and is the most 'robust' procedure. The differential version of the 
standard method relies on a stability parameter and seems to demand a prior estimate of 
the Lyapunov spectrum. The equations for tangent flow are nonlinear in this version and 
highly so in the new method. This is what makes them less efficient, though the number 
of coupled differential equations to be solved is smaller in the new method. However they 
are still useful as alternative algorithms for the computation of Lyapunov spectra. We 
have also made a comparative study of the computation of the Lyapunov eigenvectors 
using the standard method and its differential version. The eigenvectors converge fairly 
rapidly (compared to the exponents) and the two procedures yield essentially identical 
results. 
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Figure Captions 



Fig.l: Plot of the Lyapunov exponent for the for the coupled quartic oscillator system. 
The thin hne corresponds to the standard method. Of the two thick line one with '+' 
mark corresponds to the differential method and the one with the 'x' mark corresponds 
to the new method. 
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0.1738 
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(?7. = 4, a = 8) 


0.0011 


0.0011 


0.0010 


0.0001 


0.0001 


0.0001 


zi= 0.8 


-0.0011 


-0.0011 


-0.0013 


-0.0001 


-0.0001 


-0.0001 


Z2= 0.5 


-0.1892 


-0.2095 


-0.1795 


-0.1738 


-0.1793 


-0.1806 


^3= 1.0 


(0.0000) 


(0.0000) 


(0.0000) 


(0.0000) 


(0.0000) 


(0.0000) 


Z4^— 1.3 








[492.09] 


[7658.6] 


[39012.8] 



Table 1: Comparison of the Lyapunov spectrum and the computational time required 

to evaluate them with three different methods for some of the systems with n = 2, 3 and 
4. 
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Svstpm with 
initial condition 


Tjvammov f 
Standard method 


i f?pn vectors 

Differential method 


Driven van der Pol 
oscillator {n — 2) 
zi=-1.0 
Zo= 1 


Di{ 0.894, 0.447) 
D2(-0.447, 0.894) 


Di( 0.894, 0.447) 
D2(-0.447, 0.894) 


Lorenz 
system (n = 3) 
^1=0.0 
^2=1.0 

^3=0.0 


Di{ 0.004, 0.040,-0.999) 
L'2(-0. 789,-0.614,-0.028) 
D3(-0.614, 0.788, 0.029) 


Di( 0.004, 0.040,-0.999) 
D2(-0. 789,-0.614,-0.028) 
L'3(-0.614, 0.788, 0.029) 


Problem (n = 4) 
zi=1.0 
Z2=2.0 

Z3 = l.O 

Z4=0.5 


Di{ 0.230, 0.139,-0.868,-0.417) 
D2(-0.291, 0.262,-0.427, 0.815) 
D3(-0.373, 0.854, 0.187,-0.309) 
D^l 0.850, 0.427, 0.171, 0.255) 


Di( 0.233, 0.136,-0.863,-0.428) 

D2(-0.288, 0.263,-0.438, 0.810) 
D3(-0.371, 0.855, 0.188,-0.309) 
D4I 0.851, 0.425, 0.170, 0.256) 


R n*5*?l PT* 

hyperchaos (n = 4) 
Zi=-20.0 
;22= 0.0 

Z3= 0.0 

-1= 15.0 


Di{ 0.660, 0.081,-0.051, 0.745) 
L'2(-0.749, 0.115, 0.022, 0.653) 
L'3(-0.014,-0.928, 0.347, 0.137) 
L'i(-0.058,-0.345,-0.936, 0.025) 


Di( 0.660, 0.081,-0.052, 0.745) 
L'2(-0.749, 0.111, 0.022, 0.653) 
D3I 0.030, 0.991, 0.005,-0.134) 
Di( 0.050,-0.003, 0.998, 0.025) 


[n = 4, a = 6) 
zi= 0.8 
Z2= 0.5 

-23= 1.0 
Z4= 1.3 


Di{ 0.687, 0.685, 0.162, 0.182) 
D2(-0.223, 0.241, 0.672,-0.663) 
D-s{ 0.670,-0.666, 0.231,-0.232) 
i:>4(-0.170,-0.173, 0.684, 0.688) 


Di{ 0.684, 0.684, 0.178, 0.185) 
D2(-0.234, 0.241, 0.666,-0.666) 
D3I 0.666,-0.666, 0.241,-0.234) 
i:»4(-0.185,-0.178, 0.684, 0.684) 


doimled niiartir oscr 
(n = 4,q; = 8) 
zi= 0.8 
Z2^ 0.5 

-^3= 1-0 
^4= 1.3 


Di{ 0.503,-0.351, 0.581,-0.535) 
D2I 0.634, 0.744, 0.080, 0.196) 
i:>3(-0.084,-0.192, 0.638, 0.741) 
D^l 0.581,-0.536,-0.498, 0.356) 


Di{ 0.503,-0.352, 0.583,-0.533) 
D2I 0.635, 0.740, 0.085, 0.204) 
i:>3(-0.088,-0.201, 0.633, 0.743) 
D4I 0.580,-0.536,-0.502, 0.351) 



Table 2: Comparison of the Lyapunov eigenvectors computed using the differential 
and the standard methods for some systems with n = 2, 3 and 4. The eigenvectors are 
at t = 1000 for the differential method and at t = 1000, 35, 150, 170, 200, 20 respectively 
form top to bottom for the standard method. 



25 



